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ABSTRACT 

Methodologies to improve disdrometer processing, loosely based on mathematical techniques 
common to the field of particle flow and fluid mechanics, are examined and tested. The 
inclusion of advection and vertical wind field estimates appear to produce significantly improved 
results in a Lagrangian hydrometeor trajectory model, in spite of very strict assumptions of non- 
interacting hydrometeors, constant vertical air velocity, and time independent advection during 
the scan time interval. Wind field data can be extracted from each radar elevation scan by 
plotting and analyzing reflectivity contours over the disdrometer site and by collecting the radar 
radial velocity data to obtain estimates of advection. Specific regions of disdrometer spectra 
(drop size versus time) often exhibit strong gravitational sorting signatures, from which estimates 
of vertical velocity can be extracted. These independent wind field estimates become inputs and 
initial conditions to the Lagrangian trajectory simulation of falling hydrometeors. 


INTRODUCTION 

Most precipitation reaching the ground in the tropics does so in the form of raindrops and 
occasionally hail. From a very simplistic point of view, these hydrometeors can be viewed as 
non-interacting and falling at a constant terminal velocity. Even though this assumption is not 
generally accurate, because of the many physical processes such as break-up, evaporation, and 
coalescence, it is nevertheless a useful assumption for the purpose of analyzing many important 
properties of rainfall. A one-way coupled two-phase system approach for computing the 
dynamics of particles in a fluid system is standard practice when two-way coupling is not 
practical. For example, the interaction of aerosol droplets in a compressed gas stream may be, to 
first order, predicted by modeling the trajectories of individual non-interacting droplets propelled 
by a gas that is unaffected by the presence of the droplets. A similar methodology can be used to 
approach analysis of hydrometeor trajectories which may be highly influenced by the fluid 
motion of the ambient air. 
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The convention adopted in this work to describe the ambient air motion is to imply a cylindrical 
coordinate system so that wind vector at every point in space is composed of a horizontal and a 
vertical component. Keeping with the meteorological convention, the horizontal wind 
component is referred to as the advection velocity, while the vertical component is referred to as 
the updraft velocity (for a positive vertical motion), or downdraft velocity (for a negative vertical 
velocity). 

D1SDROMETER GRAVITATION SORTING SIGNATURE 

The primary motivation for this work at the start was the observation and subsequent attempt to 
model gravitational sorting of raindrops in disdrometer spectra data. Gravitational sorting 
observations in disdrometer data have been reported in the literature on numerous occasions, but 
it is not a common observation simply because of the way in which disdrometer data is typically 
stored and plotted. The gravitational sorting (GS) signature is best observed when every drop 
impact measured by the disdrometer is time tagged and then displayed as a scatter plot of drop 
diameter D vs. time /. The resulting D-t plots often show marked diagonal features, where these 
GS features are always characterized by a negative slope when plotted as D vs. t. This can be 
explained by a simplistic model of non-interacting drops traveling at terminal velocity v D as a 
function of drop diameter D, along with an ambient vertical wind motion w. At a time t = t 0 , a 

pulse of rain characterized by a typical drop size distribution (DSD), such as the Marshall Palmer 
(MP) (Marshall, 1948) or gamma distribution (see Appendix A), starts at a height h above the 
ground site where a disdrometer is acquiring raindrop spectra (i.e., counting and measuring drop 
sizes). The height h is an idealized point of rainfall generation where a rainfall DSD, such as the 
MP DSD at z - h, completely determines the size distributions of falling drops. This point 
should be correlated to a level within a cloud where rainfall begins its unimpeded trek to the 
ground, such as the 0° C isotherm. 

Other hydrometeor p can be treated similarly, but some significant difficulties with mixed phases 
(rain, hail, etc.) are encountered because of factors such as overlap in size distributions (small 
hail may be larger than large rain drops); differences in terminal velocities (round hail versus 
flattened rain drops); and differences in radar reflectivity (ice has a somewhat higher radar 
reflectivity than rain for equivalent scattering cross-sections). 

The GS disdrometer signature can be modeled based on the simplifying assumption of non- 
interacting drops by first expressing the hydrometeor fall time as: 

r D =h/(v D -w) (1) 

where r D -l D -t 0 is the fall time of a hydrometeor of diameter D, t D is the clock time of the 
disdrometer recorded hit, v D is the terminal velocity, and w is the vertical air motion. An 
estimate of terminal velocity reasonable for drop sizes associated with normal rain, is v D » aD h , 
where a « 4.5 [m s' 1 mm' b ], and b = \/2. Substituting this expression into Equation (1) and 
solving for the drop diameter D, results in: 


D = 


\ 1/6 


' hi ' t u + w 

< a , 

undefined 


h / t d + w > 0 
h/ t d + w<0 


( 2 ) 
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Figures 1 and 2 show Equation (2) plotted for various values of h and w. Equation (2) is an 
idealized description of gravitational sorting raindrops due to a zero width rain pulse (delta 
function). Figure 3 is a Monte Carlo simulation of raindrops generated by 60 [s] pulse width MP 
DSD. The simulated disdrometer data of Figure 3 differs from the real disdrometer data of 
Figure 4 for several reasons, to be discussed later. 

The vertical feature shown in the h = 0 case on the left side of Figure 3, can also occur for 
w ->- oo , where the raindrops generated at h and t = t 0 are swept downward by infinite 

velocity downdraft to appear immediately at the ground to be recorded by the disdrometer. 
Neither the h = 0 or w — > -oo cases are physical cases, but serve as a convenient means to test 
and exercise the simulations. 

Pulsed rainfall is commonly found in convective cells and thunderstorms where large variations 
of the time dependence of rainfall rate are common, where the most extreme case is a square 
pulse in both time and space. The GS signature is less common in the case of stratiform rainfall 
because the time derivative of rainfall rate is small so that no pulse-like conditions are observed. 
The D-t character shown on the left side of Figure 3 (the h = 0 case), occurs commonly in real 
disdrometer data when the rainfall rate is constant. It is also somewhat common to see the sharp 
discontinuities of D-t as shown in this simulation. This feature can occur during convective 
rainfall or stratiform rain due to advection. Advection causes the rainfall source, which is of 
finite extent, such as a convective cell, to pass over the disdrometer. In that case, the sharpness 
of the start and stop of the D-t scatter plot is due the sharpness of the spatial extent of the rainfall 
source. 

Figure 4 shows disdrometer data from the UCF Joss site in Orlando, FL (lat: 28.6016, Ion: - 
81.1 986), corresponding to September 30, 2008, from approximately 1 9:00 - 2 1 :00 UTC. Even 
thought the individual drops impacts are not time tagged in the Joss data, the stored histogram 
size interval and time interval are sufficiently small to generate a D-t scatter plot. The count in 
the histogram drop size-time interval is plotted as a density plot in Figure 4 such that a higher 
count in the histogram bin appears as a darker spot on the plot. Near the end of the rainfall 
event, GS features due to pulsed rain become clearly evident. Equation (2) is manually fitted to 
several of these features and overlaid on the plot. 

DISDROMETER DERIVED RAINFALL PRODUCTS 

From the disdrometer histogram, an equivalent drop size distribution can be immediately 
computed: 


Drop Size Distribution : 

N ik = ^ 

v D A s AtAD 

(3) 

Rainfall Rate : 

n M 

R k = S\D)H ik 

6A s Atjf J Jk 

(4) 

Reflectivity : 

1 M D 6 

Z k = f 1 H Jk 

A s Atjf lVDj Jk 

(5) 
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where H Jk is the Joss disdrometer generated histogram corresponding to elapsed time kAt and 
drop size Dj ; A s is the Joss sensor area equal to 50 [cm 2 ]; and v D is the terminal velocity (see 
Appendix) corresponding to drop diameter size Dj . In Figure 4, A/ = 1 0 [s] and the histogram 
has been resampled to a uniform diameter size so that A D = 0. 1 [mm] and D ] = jAD . In Figure 
5, Al = 30 [s] and is based on the Joss table of M= 127 non-uniformly spaced diameter sizes. 

RADAR DATA 

The Melbourne NWS radar is located approximately 55 [km] to the southeast of the UCF Joss 
site. For purposes of this study, radar data in only the lowest four scan elevations will be 
processed, and only in a 4x4 [km] horizontal extent centered over the Joss site. This region is 
later expanded to 8x8 [km] in order to extract the advection velocity. For the majority of this 
analysis, the radar reflectivity data is plotted in a pseudo 3D format, referred to as 3DRadPlot 
(see Appendix B). This display fonpat was adopted and utilized in the study of hail events 
surrounding the Pad 39B structure at the Kennedy Space Center (Lane, 2008). The same format 
is useful in analyzing the UCF Joss disdrometer data. 

The basis of 3DRadPlot display processing is based on a simplified and customized version of 
Shepard’s interpolation formula (Shepard, 1969). For any point in space, the interpolated 
reflectivity Z(r) at r = {x,y,z} is due to the weighted average of all reflectivity data points Z, at 
the center of all radar bin locations r, in the local region: 

izM- r,r 

z ( r) = -*L — (6) 

F~«i| 

where N is the number of reflectivity data points used in the interpolation, for example, N= 48, 
the number of circles in Figure 6 (the number of radar bins in the region defined by 4x4 [km] x 
4 scan elevations). The exponent parameter p controls the rate of transition between actual data 
points and the interpolated values. The larger this value, the sharper the transition between actual 
values defined by r, . In this work, a value oip = 8 was used throughout. Figures 7 and 8 are 
examples of 3DRadPlot outputs of reflectivity for 19:23 UTC and 19:38 UTC respectively. This 
can be compared to the reflectivity derived from the Joss as shown in Figure 5 (more on this 
topic below). 

JOSS DERIVED PRODUCTS 

Figure 9 is a Joss histogram corresponding to Figure 4, summed over all k time intervals. This in 
itself is not a terribly useful product, but it is easy to generate and may be serve as a convenient 
method of categorization of rainfall events. The formal drop size distribution specified by 
Equation (3) is a much more useful quantity in comparing rainfall events, especially since it is 
normalized by time. Figure 10 is the Joss derived drop size distribution using Equation (3), for 
the entire storm event of Figure 4, with At = 6000 [s]. The MP DSD equivalent corresponds to 
an average rainfall rate of 10 [mm h' 1 ]. Since the length of the storm event may be ambiguous, it 
is preferable to define shorter time intervals and then generating multiple DSDs as a function of 
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time. Of course the classic method promoted in the work of Marshall and Palmer is to group the 
DSDs by rainfall rate as opposed to time. 

Figure 1 1 is a comparison of the Joss derived rainfall rate from Equation (4) At = 30 [s] (thick 
line) and the rainfall rate from three collocated rain gauges (thin lines). The Joss and three rain 
gauges are clustered together within a 5 [m] diameter area on the roof of the UCF Engineering 
building in Orlando, FL. The total storm accumulated gauge rainfall average is 19.39 [mm], as 
compared to the 16.26 [mm] for the Joss. The cause of this 16% difference in rainfall 
accumulation between the Joss and rain gauges is not known. Suspected causes may be 
interfering wind on the roof of the building under conditions of high rainfall rate, or a 
disdrometer that is in need of calibration. This discrepancy is larger than it should be, but it is 
not the topic of this work and will not interfere significantly with the results presented in this 
paper. In future work however, we would hope to have a better initial correlation between 
disdrometer derived rainfall rate and collocated rain gauge data. 

Figure 12 is close-up of disdrometer derived reflectivity for the lowest four elevation scans over 
the UCF site (refer to Figure 6) with time marks for 19:23 UTC and 19:38 UTC scans. The 
average disdrometer Z for 19:23 UTC is 50.1 dBZ and 43.9 dBZ for 19:38 UTC when averaging 
over the two sets of four elevation scan time marks. The radar reflectivity of Figure 7 for 19:23 
UTC is ~ 45 dBZ, about 5 dBZ lower than the disdrometer derived value. The radar reflectivity 
shown in Figure 8 for 19:38 UTC is ~ 30 dBZ, about 14 dBZ lower than the disdrometer derived 
value. Proposing and testing and explanatyion for these kind of discrepancies is the goal of this 
work. The remainder of the material presented in this paper addresses this issue. The solution 
strategy takes into account the GS observations in the disdrometer data. It also should be noted 
at this time that final strategies for improving the disdrometer derived radar reflectivity will 
involve tactics that will ultimately be used to compute the disdrometer extrapolated 4D-DSD (a 
function of x, y, z, t, and D). Then, all disdrometer products can be computed from the 4D-DSD, 
including the equivalent radar reflectivity based on the 6 th moment of the 4D-DSD. 

INCORPORATION OF DROP FALL TIME 

Figure 13 depicts three methodologies for computing disdrometer derived reflectivity. The first 
method in Figure 1 3a is simply the direct application of Equation (4), as was demonstrated in the 
previous section for the 19:23 UTC radar scan, in which case a 5 dBZ discrepancy was observed 
between the disdrometer and radar reflectivites. Figure 13b illustrates an improvement to this 
simple case in which a time delay is applied to compensate for the average fall time from the 
height of the center of the radar bin based on an average drop terminal velocity. This 
intermediate case is presented without an example because our real interest is the third case of 
Figure 13c. This is the GS time delay case where the disdrometer time delay is based on 
gravitational sorting from Equation (1) and is dependent on the terminal velocity v D of each 
drop size as well as the vertical air motion w. The examples to be examined will now be based 
on comparing the simple case of Figure 13a and two case from Figure 13c, one with w = 0 and 
one with w chosen to enforce a good comparison (in a least squares sense) between the 
disdrometer Z' and the radar Z. In order to incorporate Equation (1), the computation of 
disdrometer derived reflectivity, from Equation (5), must be modified: 
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(7a) 


l M D 

z k =— — y A// 

A s& k* v D 
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k' = k + 


Int 

< 

z 

(v D; -w)At 


undefined 


v D> >w 


v Dj <w 


(7b) 


The presence of v D in Equation (7b) implies that the index k' in Equation (7a) is a function of 
index j, which depends on the details of the terminal velocity approximation used. Note the 
vertical air velocity is not used in Equation (7a) since it is assumed that the vertical air motion is 
zero near the ground, and up the height that is required to reach terminal velocity. So then for 
the largest drop sizes, it is assumed that vertical air motion is zero within a vertical region 
extending to 10 to 20 [m] above the disdrometer. Otherwise, Equation (7a) would need to be 
modified to account for the effect of vertical air movement on drop velocities at the disdrometer. 

The summation must be constrained on several fronts. To examine this in more detail, we will 
need to pick an explicit form of v D . Using the drop terminal velocity from a previous section, 

v D » aD b , and replacing D with the discrete value D J = jAD , results in: 


Int 

z 

(aAD b j b - w)At 

< 


undefined 


j > jo 

j * j 0 


Equation (7a) may now be expressed in terms of the drop size index j as: 


( 8 ) 


z k = — — ^ (yAD) H k . 
A s aAt % o (jAD) b Jk 

a M 

A S aAt Mo 


(9a) 


where. 


jo - * + 


Int 

1 

1 

a 1 * AD 


0 


w > 0 
w < 0 


(9b) 


and since the sum in Equation (9a) is truncated for j < j 0 , the time delayed index k' can be 
written as: 


k' = k + Int 


(aAD b j b - w)At 


(9c) 
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Equations (9) define an algorithm to compute a disdrometer derived reflectivity that accounts for 
drop fall time due to terminal velocity as well as vertical air motion, corresponding to the case 
depicted in Figure 13c. Equations (9) an be forced to reduce to the simple case of Figure 13a and 
Equation (5), by letting w -» -oo . The result is that k' -» k and j 0 -> 1 , so that Equation (9a) 

reduces to Equation (5). Physically, this would represents a hypothetical infinite downdraft 
which transports the drops to the ground at the disdrometer in zero time, so that the radar 
reflectivity then corresponds to the disdrometer derived reflectivity. This absurd hypothetical 
case highlights the equally absurd basis of applying the simple case of Figure 13a and Equation 
(5) for comparing disdrometer reflectivity to radar reflectivity. 

An additional limit that needs to be noted is due to the finite length of Hjk in the k direction 
(discrete time axis), k' from Equation (9c) can exceed this extent under several conditions. 

When this occurs, the summation in Equation (9c) must simply be truncated or Hjk set to zero for 
those values of k' . 

Table 1 compares the UCF Joss disdrometer derived reflectivity corresponding to consecutive 
radar scans from 19:04 - 19:38 UTC, on September 30, 2008. The 3DRadPlot format is utilized 
for easy comparison of the volume scan reflectivity in a 2x2x1 [km] box centered over the Joss 
disdrometer. Note that because of the geometry, only the lowest elevation scan will have a 
significant impact on the displayed reflectivity, even though all points in the lowest four scans 
are used (as illustrated in Figure 6) in the Shepard (1968) interpolation. Recall that the 
disdrometer data is resampled to a constant drop size bin width, with AD = 0. 1 [mm ] , and the 
original sample time, A / = 10[s] are used in Equations (9). 

In order to obtain a more meaningful comparison with radar, disdrometer Z is computed via 
Shepards interpolation at all of the radar bin points shown in Figure 6, with the actual height 
values z and time of scan defined by index k. Note that the exact positions of those points are 
contained in the Level III product data file and change position from volume scan to volume 
scan. The first column is disdrometer Z using the simple case of Figure 13a and Equation (5), 
however it is calculated by the numerically equivalent method of setting w = -999 in Equations 
(9). The second column is computed similarly as the first column, but with tv = 0. The third 
column is computed by choosing a w that gives the best (very subjective) overall match to the 
radar data of the last column. In general, there is an improvement from left to right. In some 
case, for example the 19:38 UTC scan, the first to columns are drastically off from the radar data, 
but a empirical guess at tv = 6.5 [m s' 1 ] results in a good comparison. 

Consideration of horizontal advection as well as drop vertical fall time provides additional 
improvement in the calculation of disdrometer reflectivity, as will be shown in the following 
section. 

INCORPORATION OF CLOUD ADVECTION EFFECTS 

In this approach, radar reflectivity and radial velocity data are used in addition to disdrometer 
spectra, where the end goal is no longer just a improved comparison between radar Z and 
disdrometer Z. The goal is to arrive at a good prediction of the DSD as a function of all three 
spatial coordinates, as well as time. Once that is accomplished, the resulting 4D-DSD can then 
be used to calculate rainfall products, analogous to Equations A. 18 - A.26 of Appendix A. 
Never-the-less, it is still convenient and useful to compare the final disdrometer derived 
reflectivity from the 6 th moment of the 4D-DSD to the corresponding radar reflectivity. 
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Even though there are typically four dozen radar bin locations at each complete scan in the 
2x2x1 [km] volume surrounding the UCF disdrometer site as shown in Figure 6, far fewer than 
that number contribute significantly to the estimate of the 4D-DSD, primarily those in the lowest 
elevation scan. Based on Figure 6, there are a half dozen especially influential points. 

Therefore, the number of free and arbitrary model fitting parameters should not be allowed to 
exceed the number of radar reflectivity points in order to properly satisfy least squares fitting 
requirements. However, this is more of a intelligent guideline than a strict rule since the 
approach taken in this work is not to implement a stringent least squares comparison of radar and 
disdrometer Z. 

Table 2 shows the Melbourne NEXRAD radial velocity over the disdrometer site for several 
scans of the example rainfall event for the first three elevations. This data was acquired using 
the NOAA Weather and Climate Toolkit, version 2.2, from the National Climatic Data Center (a 
free software download). Note that both classified (velocity quantized to multiple of 5 [kts]) and 
unclassified data are shown in Table 2. The important characteristic to note in this data is the 
sometimes strong vertical gradient of the horizontal winds. This feature is most notable at the 
beginning of the storm where many atmospheric processes may interact, such as gust fronts and 
sea breeze collisions, as well as strong updrafts and downdrafts. On the backside of the storm 
(last two rows of Table 2), the vertical gradient is small or non-existent. 

An empirical advection model that is both relevant to the 4D-DSD volume and consistent with 
the observations of vertical advection gradients in the NEXRAD radial velocity, was used in this 
work: 


with. 


u(z)= r + (l-/)tanh^y^ju 0 


u<> 


foSlf/ 

f 0 


(10a) 


(10b) 


The model described Equation (10) exhibits two asymptotic values: u(z) — > (2y - l)u 0 for 
(z-z 0 )/ L. «0 and u(z)->u 0 for (z-z 0 )/ L. »0. The center of the transition at z = z 0 is 
u(z) = y u 0 and the rate of the transition between asymptotic values is controlled by the 
parameter L. . The 4D-DSD model is implemented in a Cartesian coordinate system with 
positive jc pointing east, positive v pointing north, and positive z pointing up. Therefore, 
advection angle if/ adheres to the polar angle convention of the Cartesian system so that if/ = 0 
simulates an advection travelling from west to east, if/ = 90° simulates an advection travelling 
from south to north, and so on. 


A simulated radial velocity can be calculated from Equation (10a) by performing the vector dot 
product with the radar direction vector: 


U r (z) 


fiosOcosf 
cos Q sin </> 

< sin# J 


Mz)^ 


w , 


(11a) 


2/27/2009 


Page 8 


s 
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where <j> direction from the radar to the disdrometer site (in the Cartesian coordinate system of 
the 4D-DSD model); 9 is the local radar scan elevation angle above the disdrometer; 9 0 is the 
radar scan elevation angle relative to the radar site; 5 is the distance along the earth’s surface 
from radar site to disdrometer site; and R E is the average earth radius, 637 1 [km]. The last term 
in Equation (1 lb) takes into account the curvature of the earth, while the factor of 4/3 accounts 
for radar refraction. Equation (lib) assumes that the radar and disdrometer sites are both at 
points on a perfect spherical earth of radius R E . 


Since under most conditions, the vertical component of U r (z ) is very small compared to the 
horizontal component, setting 9 = 0 simplifies the simulated radial velocity. Combining 
Equations (10) and (11) with 9 = 0 yields: 


U r (z)»u r (z) = u 0 \ y+(l-y) tanh 


z-z n 


cos 


( 12 ) 


In order to compare the model radial velocities from Equation (12) to the NEXRAD radial 
velocities V r , such as those in Table 2, a vertical integration needs to be performed over the radar 
beam width at each elevation angle: 


u r {m) = - 


Z m2 Z m\ z, 


f 


= u. 


L. 


/ + (1 - y) — — — logl cosh Z ° Z " 2 • sech *° " ml 


Z n ~Z„ 


\\ 


(13) 


Z ml Z m\ 


L. 


COS (0-t//) 


JJ 


where z ml and z m2 are the bottom and top edges of the beam at the mth radar beam elevation 
angle. 

The average advection velocity u(/n) the wth radar beam elevation angle can be estimated by 
plotting equal dBZ contours between consecutive scans, such that the contour lines bracket the 
disdrometer site. It is necessary to use the appropriate contour, corresponding to the dBZ value 
passing over the ground site at each scan time. For example. Figure 14 shows the NEXRAD 5 
dBZ reflectivity contour for the m = 1 scan (lowest elevation scan) centered over the disdrometer 
site for the 19:04 and 19:08 UTC scans. The procedure depicted in Figure 14 can be repeated for 
all elevation angles of interest and for all time scans of interest. 

Once the Z-contour velocities u c (m) have been estimated from the radar Z data, for at least m = 
1, self consistency of the radar data can be checked by comparing the radial components of 
u c (/w) with the Doppler radial velocities from Table 2. The radial component of u c (w)are 

determined by taking the dot product with the direction vector from the radar, similar to Equation 
(11). This comparison provides some confidence of the integrity of the radar derived advection 
data. In some cases, the Doppler radial velocity may be unreliable because and spatial and 
temporal fluctuations. The Z-contour derived advection may be unreliable in some cases, such as 
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during collisions of the sea breeze and front boundary movement, in which case advection may 
be undefined. For the purpose of disdrometer spatial and temporal interpolation-extrapolation of 
drop size distributions, the u c (/w) may be a better measure of relevant advection. For example, 

during times a convective cell is stationary, u c (/w) would most likely be a better indicator than 
the associated Doppler velocities. 

The goal of the discussion of the last several paragraphs is to choose the best parameters for the 
advection model. Equations (10), based on empirical fitting to the Doppler and reflectivity 
derived horizontal air motion at the first two or three elevations angles. For example, the 
advection model corresponding to the 19:04 UTC radar scan is shown in Figure 15. It should be 
noted that the advection model of Equation (10a) only has the ability to change the velocity 
magnitude, not the direction. It can however, reverse the sign. Increasing the capabilities of the 
advection model would be a candidate for future work. 

GENERALIZED LAGRANGIAN RAINFALL MODEL 

The previous two sections described independent horizontal and vertical trajectory models of 
falling rain drops due to gravity and the effects of ambient air motion. The primary effect of 
atmospheric air motion is due to drag forces that either move the particle along with the ambient 
air and/or limit the drop fall velocity, i.e. terminal velocity. Other effects include particle lift due 
to the rotation of the particle or rotation of the gas around the particle. Lift forces, which may be 
positive or negative, are always in the direction of the gravity vector and are proportional to the 
horizontal component of drag. Lift is usually much smaller in magnitude than drag forces, 
except in the case of very fast rotations. Mass loss or gain is another affect that controls the 
detailed motion of rain drops. Evaporation, collisions, spontaneous breakup, and coalescence are 
mechanisms that further complicate drop dynamics. 

A traditional approach to simulating trajectories of single particles in a gas flow is to employ a 
recursive integration algorithm such as the fourth order Runge-Kutta (RK4) algorithm, to 
compute the particle position at each time step. The inputs to the integration algorithm include 
local gas velocity vector, density, viscosity, and temperature at each of the particle’s time- 
stepped positions. For some applications such as aeronautics, these gas properties might be 
output from a computation fluid dynamics (CFD) or direct simulation Monte Carlo (DSMC) 
model of the gas flow. Advanced fluid mechanics models utilize a two-way coupled system 
where the gas flow is affected by the particle flow, and the particle flow is in turn affected by the 
gas flow. Less computationally intensive software systems implement a one-way coupled model 
where the gas flow affects the particle flow, but the particle flow has no affect on the gas flow. 

The horizontal advection and vertical wind components of falling rain can be combined in a one- 
way coupled algorithm, using physical property formulas of the surrounding air as input to the 
Lagrangian trajectory integration. For a hydrometeor of diameter D and mass m, the trajectory is 
due to the external forces: gravity, vertical updrafts/downdrafts, and advection. The sum of 
external forces on a hydrometeor is equal to its acceleration, which can be estimated by a Taylor 
series expansion about time point n, resulting in a set of difference equations for position and 
velocity (see Appendix C for a RK4 implementation): 

=v„_i+a„_,A/ (14a) 
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where g E is earth gravity (9.80665 [m s' 2 ]) ; p H is hydrometeor density (997.0479 [kg m' 3 ] at 
25 C); p{y k ) is the air density; and U(r Jt ) is the total air velocity at the Mi time step particle 
position r k . The direction of the gravity unit vector e £ is given by: 


e £ - 


-1 


Jx 2 +y 2 +(R e +x) 2 


x ' 


y 

R e+z, 


^o> 

o 


(15) 


for jc,y,z « R e 


The coefficient of drag, C D is a function of the Reynolds number. Re : 

Dp,\V. ~ v„| 


Re. =■ 


P« 


(16) 


where p„ is the dynamic viscosity of air and is related to temperature by: 


T + C 

ft =£(<•,)=«, ° 


T 

\ l o y 


T(r k ) + C 

where p 0 = 1.827 xlO' 5 [Pas], T 0 =291.15 [K], and C 0 s 120[K], 

The air density p„ in Equation (14d) and (16) is related to temperature and air pressure: 

Pjr m )M 0 
R'T( r.) 


(17) 


(18a) 


Page 

11 


2/27/2009 



(18b) 


T(r n ) = T(0)-rjz n 

where M 0 is the molecular weight of dry air, equal to 0.0289644 [kg mol' 1 ]; R‘ is the gas 


(18c) 


constant, equal to 8.31432 [J mol' 1 K' 1 ]; P 0 is standard air pressure at sea level, equal to 101325 
[Pa]; rj is the average rate of temperature change with altitude, equal to 0.0065 [K m' 1 ]; and 7(0) 
is the standard temperature at sea level, equal to 288. 1 5 [K], 

The coefficient of drag is computed from the following empirical formula: 


where the drop terminal velocity v D is given by an approximation that accounts for of altitude 
effects, such as Best (1950) (see Equation A. 16) in Appendix A. 

Equations (14) through (20) define a complete Lagrangian trajectory model for hydrometeors, 
one that provides a method to integrate the equations of motion, following the hydrometeor path 
from an arbitrary stating point r 0 = i.x 0 ,y 0 ,z 0 ) t0 the ground. The time of travel r = L At is 

found by noting the number of time steps L where r, — > (x, , y L ,0) . The general strategy is to 
choose a set of starting points r 0/ , then run the trajectory Equations ( 14) - (20) recursively until 
the hydrometeor hits the ground., repeating this for a all discrete values of drop diameter Dj of 
the relevant drop size range. For every D J starting at r 0i , the coordinates on the ground and the 

arrival time at the ground is logged. The times and locations of the hydrometeor ground strikes 
are then related to the disdrometer spectra by a time delay and spatial extrapolation from the 
disdrometer location to the location of the hydrometeor strike. This procedure provides a method 
to generate a volume drop size distribution (DSD) in 3D space and in time. The Lagrangian 
approach to computing a 4D-DSD from the disdrometer data is a powerful method since it can 
easily accommodates complex spatial and temporal advection functions, as well as complex 
vertical wind profile functions. In previous work, a 3D-DSD algorithm was implemented in 
FORTRAN 95 to process disdrometer data using a set of ballistic equations which had been 
integrated analytically using simple functions of advection and vertical wind motion (Lane, 
2000), where the wind motion was assumed to be time independent over the period of drop 
trajectory calculations. The Lagrangian trajectory method has no requirement of time 


24.0 Re ' Re <2 

C D =< 18.5/te” 06 2 < Re <500 
0.44 Re > 500 


(19) 


The initial conditions assume that the hydrometeor is at terminal fall velocity: 



( 20 ) 
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independent wind motion. (Note that the difference in terminology adopted in this work is that 
previous work computed a 3D-DSD where time was involved only in the trajectory equations, 
but was not used to compute changes in acceleration as a function of time). 


D1SDROMETER DERIVED REFLECTIVITY 

A final problem to solve is that of spatial extrapolation from the disdrometer site to the point on 
the ground of the hydrometeor strike. An ad hoc approach is to assume that each yth-£th bin of 
the disdrometer histogram H Jk travels from the disdrometer site to the hydrometeor impact point 

p v = ( Wl> ) using the following transformation (Lane, 2000): 


4 ~ t i + X ij ~ 


»(°)P, 

|u(0)| 2 


( 21 ) 


where u(0) is the advection velocity near the ground at z * 0 ; p,. is the point of impact on the 
ground for the yth drop size falling from the r 0i position in space; r y is the fall time of the yth 
drop size from the point r 0 . ; and /, is the absolute time that drops begin to fall from the point 
V The disdrometer data corresponding to H jk . is then used to estimate the extrapolated 
disdrometer spectra at the hydrometeor impact point p v , where: 


k' 



( 22 ) 


and t 0 is the absolute start time of the disdrometer spectra acquisition that produces H jk . . 

For the purpose of comparing disdrometer reflectivity to radar reflectivity, it is convenient to 
define the set of r 0 . and t, corresponding to the set of radar bins over the disdrometer volume, as 
shown in Figure 6. 

In some cases, such as times when a sea breeze collides with a convective cell front, advection 
may not be well defined. In such cases where advection goes to zero or reverses direction, the 
last term on the right (the advection term) in Equation (22) can be deleted, which essentially 
reduces to the case of Equation (7b). In cases where advection is non-zero, but small or simply 
not consistent, a linear combination of these two solutions can then be used to compute the drop 
size distribution: 


N(D) = p N a (D) + (l - P) N s (D) (23) 

where J3 is an empirical parameter that mixes some percentage of the advective case N a (D) 
defined by Equation (22) and the stationary non-advective case, N S (D ) of Equation (7b). 

DISCUSSION AND SUMMARY 

Referring to the plots in Tables 1 and 3, the 19:04 and 19:08 scans show a distinct improvement 
in the disdrometer derived reflectivity when vertical fall time is incorporated into the processing 
algorithm (second and third columns compared to the first column). The Lagrangian-advection 
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model, shown in the middle column of Table 3 also shows good agreement with the radar in 
most cases. Scans 19: 13, 19: 18, and 19:23 all show good agreement in all columns with the 
corresponding radar reflectivity. One explanation for this behavior is that during the 19:13 — 
19:23 period, rainfall is characterized by a very broad pulse, so that gravitational sorting is not a 
significant effect during that time interval (refer to the disdrometer spectra plot of Figure 4 at 780 
- 1380 [s]). In this case, the disdrometer reflectivity can easily track the radar reflectivity. 

Scans 19:28, 19:33, and 19:38 begin to show severe disagreement between disdrometer and 
radar. Only the processing methods shown in Table 4 display reasonable agreement with the 
radar reflectivity. An explanation for this case is that the rainfall is characterized by impulsive 
conditions (refer to Figure 4 at 1680 - 2280 [s]). During this time, advection begins to 
degenerate due to interaction of the storm movement with surrounding opposing winds (this is 
somewhat apparent in the radar data). Another symptom of poorly defined advection is seen in 
Table 4, by examining the evolution of the mixing parameter (3. Recall from Equation (23) that 
p = 1 is a purely advective solution, whereas p = 0 is a purely non advective solution defined by 
Equation (7b). In Table 4, where all parameters are created by a manual process of finding a best 
fit to the corresponding radar plot, (5 = 1 during the non-impulsive phase, then begins to 
decreases to smaller values. The only exception, p during the 19:04 scan is less than 1, which 
may be explained by poor advection conditions at the earliest approach of the storm, when 
rainfall rate and reflectivity are very small. 

Another interesting feature in the parameter set of Table 4 is the evolution of the vertical velocity 
parameter w. Vertical velocity starts out at a small value or negative (downdraft), then begins to 
rise to larger and larger positive values (updraft), where the maximum coincides with the 
maximum radar reflector at 19:23. 

Methods to improve disdrometer processing, loosely based on mathematical techniques common 
in the field of particle flow and fluid mechanics, have been explored. The inclusion of advection 
and vertical winds appears to produce significantly improved results, in spite of very strict 
assumption of non-interacting hydrometeors, constant vertical air velocity, and time independent 
advection during the scan time interval. Future work should focus on relaxing these assumptions 
by developing at the least a simple model of drop interaction. Time dependent advection may be 
incorporated quite simply within the frame work of the Lagrangian trajectory mechanics. 

A most important activity of future work should be to exercise the model vigorously with a 
sufficient volume of data so that statistics of the model performance can be quantified. Along 
with this effort, all available independent wind field data should be rigorously collected and 
processed to be used as inputs into the Lagrangian model. Sources of this data may be derived 
directly from each radar elevation scan by plotting and analyzing reflectivity contours over the 
disdrometer site, such as that demonstrated in Figure 14, and by collecting the radar radial 
velocity data. Strong GS signatures in the disdrometer spectra are another potential source of 
vertical wind data. Other sources of data that would be very useful, if available - nearby 
anemometer or wind tower data. Rain gauge data from multiple gauges positioned within a 
kilometer or less of the disdrometer site would provide a valuable addition to the set of data 
required 4D-DSD optimization and verification. 
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0 500 1000 1500 2000 

t-tO [si 

Figure 1 . Plots of D(l) due to idealized pulse rain characterized by a delta function, using 
Equation (2) for various values of h with w = 0. 




Figure 2. Plots of D(t) due to idealized pulse rain characterized by a delta function, using 
Equation (2) for various values of w with h = 2000 [m]. 
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Figure 3. Simulated disdrometer data for a r= 60 [s] pulse of R = 100 [mm/h] rainfall using 10 6 
Monte Carlo drops generated with a MP drop size distribution. 
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1000 2000 3000 4000 5000 
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Figure 4. Joss disdrometer data from the UCF site, September 30, 2008, with / 0 = 19:00 UTC. 
D(t) plotted is from Equation (2): dashed line, w = -1.5 [m/s], h = 4000 [m]; thick dashed line, w 
= -1.5 [m/s], h = 4000 [m]; dot-dashed line, w = -1.0 [m/s], h = 4000 [m]; solid line, w = +0.5 

[m/s], h = 2000 [m]. 
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Absolute Time [UTC] 

Figure 5. Rainfall products derived from September 30, 2008 Joss disdrometer at UCF: thin line 
is rainfall rate computed by Equation (4); thick line is equivalent radar reflectivity computed 

using Equation (5). 



Figure 6. UCF disdrometer site (Joss disdrometer in located in center). Circles represent 
locations of the Melbourne NWS radar bins for the lowest four elevations (09-30-08): black 
circles, z = 1 .00 [km], / = 19:24:06 UTC; brown circles, z = 2.40 [km], t = 19:24:41 UTC; red 
circles, z = 3.50 [km], I = 19:25:23 UTC; and yellow circles, z = 4.90 [km], / = 19:25:45 UTC. 
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Figure 7. 3DRadPlot of Melbourne radar reflectivity over UCF Joss site, 09-30-08, 19:23 UTC. 




Figure 8. 3DRadPlot of Melbourne radar reflectivity over UCF Joss site, 09-30-08, 19:38 UTC. 
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1400 



Figure 9. Joss histogram corresponding to Figure 4, summed over all k time intervals. 



D [mm] 

Figure 10. Joss derived drop size distribution using Equation (3), for entire storm event of 
Figure 4, with A / = 6000 [s]. The MP DSD equivalent corresponds to an average rainfall rate of 

10 [mm h 1 ] 
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t [hr] 

Figure 1 1. Thick line is rainfall rate derived from the Joss data using Equation (4) with 
A t = 30 [s]. The three thin lines are rainfall rate measured by three collocated rain gauges. The 
total storm accumulated gauge rainfall average is 19.39 [mm], as compared to the 16.26 [mm] 

for the Joss. 



Absolute Time [UTC] 

Figure 12. Close-up of disdrometer derived reflectivity corresponding to the lowest four 
elevation scans over the UCF site (refer to Figure 6) for 19:23 UTC and 19:38 UTC scans with 
time marks. The average Z for 19:23 UTC is 50.1 dBZ and 43.9 dBZ for 19:38 UTC. This 
should be compared to the reflectivity values shown by the NEXRAD data of Figures 7 and 8. 
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(b) 


t' = t + h/ v n 


* 



(c) 


/' = / + h/(v n - vi’) 


Figure 13. Comparison strategies for radar reflectivity Z and disdrometer derived reflectivity 
Z' : (a) simple case where the time / of the radar scan is matched to the corresponding 
disdrometer data at (' = t ; (b) the simple time delay case where Z(/) is matched to the 
disdrometer data delayed by the average drop terminal velocity v D falling from a height z = h ; 
(c) GS lime delay case where the time delay is dependent on the terminal velocity v D of each 

drop size as well as vertical air velocity w. 
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Table 1. Application of Equations (9a)-(9c) to compute disdrometer derived Z with comparison 
to Melbourne equivalent radar reflectivity on a scan by scan basis, using UCF Joss data of 
September 30, 2008. Note that the simple case shown in the first column was actually produced 

by setting w = -999 [ms' 1 ] in Equations (9). 


Disdrometer derived 
reflectivity Z using 
simple case of Figure 
1 3a and Equation (5). 


Disdrometer Z using 
method of Figure 13c 
and Equations (9) 
with w set 0. 


Disdrometer Z using 
method of Figure 13c 
and Equations (9) 
with manual 
optimization of w. 


Melbourne NEXRAD 
reflectivity plotted in 
a 2x2x1 [km] volume 
centered over UCF 
Joss site. 
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w = 0 



w = +4.5 [m s' ] 



19:23 UTC 



19:28 UTC 




w = 0 


w = -1 .0 [m s' 1 ] 



19:28 UTC 


- i - 






19:33 UTC 


« 




w = 0 







w = -1.5 [m s' ] 




19:38 UTC 



w = 0 



w = 6.5 [m s' ] 


19:38 UTC 
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Table 2. Level II Melbourne NEXRAD radial velocity data over 
the UCF Joss disdrometer site, Septmeber 30, 2008. 


/ UTC 

Elev 

[deg] 

z[m] 

Classified 
V r [m/s] 

Unclassified 
V r [m/s] 


0.47 

959 

- 

-2.52056 

19:04 

1.44 

2246 

- 

-1.49176 


2.39 

3511 

-10.288 

-5.50408 


0.47 

960 

-5.144 

-4.47528 

19:08 

1.44 

2246 

5.144 

1.49176 


2.39 

3513 

-5.144 

-1.49176 


0.47 

960 

-10.288 

-6.01848 

19:18 

1.44 

2247 

5.144 

0.97736 


2.39 

3512 

9.7736 

2.52056 


0.47 

962 

-10.288 

-8.02464 

19:23 

1.44 

2245 

-5.144 

-2.52056 


2.39 

3515 

5.144 

0.5144 


0.47 

963 

-10.288 

-6.99584 

19:28 

1.44 

2246 

5.144 

1.49176 


2.39 

3508 

5.144 

0.97736 


0.47 

963 

-5.144 

-4.47528 

20:22 

1.43 

2243 

-5.144 

-4.98968 


2.39 

3513 

-5.144 

-3.49792 


0.47 

959 

-5.144 

-4.01232 

20:27 

1.43 

2242 

-5.144 

-4.47528 


2.39 

3515 

-5.144 

-4.01232 
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E-W Distance [km] 

Figure 14. Lowest elevation angle plot of the Melbourne NEXRAD 5 dBZ reflectivity contour 
over the UCF disdrometer site. Based on the translation of the equal contour line, the 
corresponding advection velocity over the site, using the notation of Equation (13), is 
approximately w c (l)« 3.3 [m s' 1 ], with y/= -68°. 




Wind Velocity 


Figure 15. Advection model for scan 19:04 UTC from Equations ( 10), with ^ = 1.4, 
/,.= 700[m], z 0 = 1600[m], y/ = -64° , u 0 = 1.6 [m s' 1 ], with <f> = 134.5°. 
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Table 3. Comparison of Lagrangian trajectory model disdrometer derived Z with Melbourne 
equivalent radar reflectivity on a scan by scan basis, using UCF Joss data of September 30, 2008. 
The first column is copied from the third column of Table 1. 
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w = +4.5 [m s' 1 ] 




■ 

w = -1.5 [m s' 1 ] 
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Table 4. Lagrangian model parameters used to produce second 
column in Table 3. 


UTC Time 

19:04 

19:08 

19:13 

19:18 

19:23 

19:28 

19:33 

« 0 [m s' 1 ] 

1.6 

-1.3 

6.2 

-1.8 

1.35 

-1.6 

5.5 

¥ [deg] 

-64 

-99 

-72 

-43 

-69 

-35 

^to 

fi 

0.61 

1.0 

1.0 

1.0 

1.0 

0.55 

0.29 

7 

1.4 

-4.0 

0.52 

-7.7 

4.2 

-7.0 

1.0 

Mm] 

700 

920 

350 

800 

850 

600 

500 

z 0 [m] 

1600 

0 

0 

0 

1500 

600 

0 

w [m s' 1 ] 

0 

0 

-2.5 

1.0 

4.5 

4.0 

-5.7 
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APPENDIX A 


DROP SIZE DISTRIBUTION 

The drop size distribution is the number of raindrop hydrometeors per unit volume per drop size 
diameter D. Usually the DSD is denoted by N(D) and is expressed in units of m' 3 mm’ 1 or cm -4 . 
The total number of drops N T in a volume of air is equal to the integral of N(D) over all drop 
diameters: 


N T = \N(D)dD (A. I) 

0 


where N(D)dD is the number of drops per unit volume with diameters between D and D + dD. 

Since the largest raindrop sizes found in nature do not exceed 6 to 7 mm, the DSD is often 
defined as: 


\N(D) 


D<D 


MAX 


DSD = \ 


(A.2) 


0 


D> D 


MAX 


Quantities involving the DSD, such as Nt , are then found from the truncated integral: 

^ MAX °° 


N t = \N(D)dD*\N(D)dD 


(A.3) 


For most quantities involving integrals of N(D), the error associated with the choice of large 
drop size cutoff is insignificant as compared to errors from other sources. A practical strategy, 
which will be used in this work, is to limit large drop sizes to 6 mm when performing data 
analysis. When performing analytical calculations on integrals involving A \D), it is more 
convenient to use infinity as the upper bound for drop size since in this case, the integrals often 
reduce to the form of a simple gamma Junction. 


EXPONENTIAL DISTRIBUTION 

The simplest and probably most commonly used DSD is the exponential distribution : 

N(D)=N 0 e~ AD (A.4) 

where N 0 and A are fitting parameters. In this case, the total number of drops per unit volume 
reduces to: 

oo 

N t =\N 0 e™ dD 

o (A.5) 

= No A' 1 
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GAMMA DISTRIBUTION 


The gamma distribution (Ulbrich, 1983) was proposed in order to alleviate a problem at small 
drop diameters seen by researchers, which suggested that N(D) — » 0 as D — > 0 : 

N(D) = N 0 D v e~ AD (A.6) 


The gamma distribution also allowed for a three parameter fit, where v , as well as N 0 and 
A are the fitting parameters. In this case, the total number of drops per unit volume reduces to: 


N t = J N 0 D V e _AD dD 

0 

= r(l + v) 

A V ) 


(A.7) 


Comparing Equations (A.4) and (A.5) to Equations (A.6) and (A.7), it can be seen that the 
exponential distribution is a special case of the gamma distribution where v = 0. 


MARSHALL-PALMER DISTRIBUTION 


The Marshall-Palmer (MP) DSD is a special case of the exponential distribution, described by 
Equation (A.4), where historically the exponential DSD may have resulted as a generalization of 
the MP DSD. The MP DSD is characterized by its explicit definition of the fitting parameters, 
N 0 and A: 


N 0 = 8000 m'^mm' 1 
A = 4.1 mm' 1 


(A.8) 


with R in units of mm h' 1 . The original DSD data by Marshall and Palmer is shown in Figure 
A-l. 

The total number of drops for the MP DSD, based on Equation (A.5), is: 

N t = N 0 A -1 

= (N 0 /4.\)R 021 (A 

= (1951.2) R 021 m* 3 

According to the MP DSD and Equation (A.9), the total number of drops per unit volume 
increase as the 0.2 1 power of R. 


RAINDROP TERMINAL VELOCITY 

The equation of motion for a raindrop of mass m D in a uniform gravitational field of 
acceleration g, falling in the z direction (defined positive from cloud to ground, where z = 0 at 
the initial starting point in the cloud), is: 
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m D z = m D g-p D z 2 


(A. 10) 


where the first term on the right-hand-side (RHS) is the gravitational force and the second term 
on the RHS is the force due to frictional drag through air. Approximations for the quantities in 
Equation (A. 10) can be estimated by assuming a spherical drop shape of diameter D and by 
neglecting the change in atmospheric density with altitude: 


m D 


6 


(A. 11 a) 


M D 


8 


(A. lib) 


where pw is the density of water, p is the density of air (neglecting altitude dependence), and 
Cd is the drag coefficient of a spherical drop. By substituting v = z in the equation of motion, 
Equation (A. 10) becomes: 


dv 
v — 
dz 


g- — V 2 

m D 

g(l-v 2 /v£) 


(A. 12) 


where v D s ^]m D g / p D . The solution to the differential equation of Equation (A. 12) is: 


v = 


v D (l 


(A 13) 


As z — ^ oo in Equation (A. 13), v — > v D , thus implying that v D is the terminal velocity of a 
drop of diameter D. Experimental data (Laws, 1941) for drop velocities versus fall height for 
several drop sizes, is shown in Figure A-2, along with the velocities predicted by Equation 
(A. 13). Note that Equation (A. 13) also correctly predicts that as z goes to zero: 

V ( 2 £ Z )' /2 (A. 14) 

which is the free fall solution in a vacuum, where the drag force vanishes. 

One of the primary problems with the former arguments is the assumption that the raindrop 
shape is spherical. Figures A-3 and A-4 show high-speed photographs of drops at terminal 
velocity (Edgarton, 1939, and Pruppacher, 1970). For drop diameters larger than about 3 mm, 
the shape becomes more flattened as the diameter increases. For this reason, terminal velocity 
formulas are usually empirical and are not usually based on the spherical drop shape assumption 
of Equation (A. I lb). 


APPROXIMATION BY GUNN (1948) 

A simple form of drop terminal velocity v D is described by a simple power-law (Gunn, 1948): 

v d *K g D V2 (A. 15) 
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where Kg = 4.5 when drop diameter D is expressed in mm and terminal velocity is in units of 
m s'*. This particular form of v D is particularly useful when evaluating integrals over drop size, 
involving drop size distribution functions. Specific examples of this will be presented in sections 
to follow. 


APPROXIMATION BY BEST (1950) 

An improved approximation for terminal velocity (Best, 1950), which takes into account the 
variation of velocity with altitude, due to atmospheric density, is: 

v D l-e (0/<,) "} (A. 16) 

where A = 9.58 m s' 1 , b = 0.0854 km' 1 , z is the elevation at the ground (or at the point of 
measurement), D is the drop diameter measured in mm, a = 1.77 mm, n = 1.147, and v D is 
expressed in units of m s' 1 . 

APPROXIMATION BY ATLAS (1973) 

Another approximation for terminal velocity (Atlas, 1973), which is also an improvement over 
the power-law form of Gunn, is: 

v D *v 0 (l -Ke cD ) (A.17) 

where v 0 = 9.65ms' 1 , K = 1.067, c = 0.6 mm' 1 , and D is expressed in mm. Even though this 
form of v D is not as simple as Gunn’s approximation, it nevertheless results in a form which is 
relatively easy to manage when evaluating integrals involving DSD functions. 


MOMENTS OF THE DSD 

Many useful atmospheric and meteorological quantities can be calculated simply by evaluating 
the xth moment of the DSD: 


M,=\D*N(D)dD (A . 1 8) 

0 

For example, the total number of drops Nt in a volume of air is equal to A/ 0 , the zeroth moment 
of the DSD, as verified by Equation (A.l). 


LIQUID WATER CONTENT 

The third moment A /3 of the DSD is proportional to the amount of water in a unit volume of air: 


Page 

33 


2/27/2009 


(A. 19) 


V = fp r \D 3 N(D)dD 


0 


= iPwM, 

If the exponential DSD described by Equation (A.5) is used to evaluate Mi , the result is: 


where pn- is the density of water. 


RAINFALL RATE 

The integral over D of the product of terminal velocity and DSD is equal to the rainfall rate. If 
the terminal velocity v D is described by Gunn’s approximation. Equation (A. 15), rainfall rate R 
becomes the 7/2 moment of the DSD: 


oo 


e' w dD 


0 


(A.20) 


— K PiyN q A 



(A.21) 


o 



where K G = 4.5 (drop diameter D is expressed in mm and terminal velocity in units of m s' 1 ). 
If the exponential DSD of Equation (A.5) is used, R then evaluates to: 


R = (0.0036)f K g N 0 j D 7,2 e~ AD dD 


(A.22) 


o 


= (0.0036)§ tt 3,2 K g N 0 A~ 9 ' 2 

where the scale factor 0.0036 places R in standard unit of mm h' 1 . 


RADAR REFLECTIVITY 


The sixth moment M6 of the DSD is approximately equal to the radar reflectivity: 


oo 



(A.23) 


o 


If the exponential DSD of Equation (A.5) is used, Z then evaluates to: 
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(A.24) 


GO 

Z=N 0 \D 6 e~ AD dD 

0 

= 720 N 0 A~ 7 

where Z is in standard units of mm 6 m' 3 . Note that a commonly used measure of Z is dBZ = 10 
logio(Z). 


Z-R RELATION 


The rainfall rate R described by Equation (A.22) can be combined with the radar reflectivity Z 
given by Equation (A.24), by eliminating A: 


Z = 


r 274211 "I 

*-14/9 x r5/9 

Va g N 0 ) 


R W9 


(A.25) 


Using K g = 4.5 and N 0 = 8000, Z becomes: 

Z = aR b (A.26) 

where a « 179.3 and 6 = 14/9 « 1.56 . Equation (A.26) is a Z-R relation derived from Gunn’s 
terminal velocity and the exponential DSD. Equation (A.26) is the well known Z-R power law 
used to convert weather surveillance radar (WSR) data to estimations of rainfall amounts. 
However, it should be noted that the National Weather Service (NWS) has adopted the modified 
values of a = 300 and b = 1.4 for use with its WSR-88D radar network. 
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APPENDIX B 

3D-MAPPING OF NWS REFLECTIVITY DATA 


3DRadPlot was originally developed to compare NWS NEXRAD volume reflectivity data to 
drop size distribution data acquired from a cluster of disdrometers. In this work, particular 
attention was given to studying data from a cluster of hail disdrometers and the corresponding 
radar data: http://www.nasa.gov/mission pages/shuttle/behindscenes/hail monitor.html . 
Additional information concerning spatial interpolation and extrapolation of disdrometers is 
included in the \Papers directory from the 3DRadPlot page at openchannelsoftware.com. 
However, 3DRadPlot plots radar data only, even though the development was motivated by 
plotting the 3D data from the disdrometer cluster. 

NCDC NEXRAD DATA INVENTORY 

Radar data from all NWS NEXRAD sites is archived at the National Climatic Data Center 
(NCDC). That data can be readily accessed at http://www.ncdc.noaa.gov/nexradinv/ . 

3DRadPlot plots Level III reflectivity data at four scan elevations. These files are denoted by the 
sub-string, NOR, N 1 R, N2R, and N3R. A complete file name is for example: 

7000KMLB_SDUS52_N0RMLB_200702262214 . txt 

corresponding to Melbourne radar data, lowest elevation scan, of February 26, 2007, 2214:00 
UTC. The additional three elevation scans that are available from Level III data, for this 
particular case, turn out to be: 

7000KMLB_SDUS22_N1RMLB_20070226221 4 . txt 
7000KMLBSDUS22N2RMLB20070226221 4 . txt 
7000KMLB_SDUS32_N3RMLB_200702262214 . txt 

These are binary files. After ordering the data from the date and time of interest, the following 
message will indicate that data is ready for pickup. The files are then downloaded by clicking on 
the link in an email message, such as the following: 


From: Orders HAS [mailto:orders@thunder.ncdc.noaa.gov] 

Sent: Sunday, August 26, 2007 2:28 PM 
To: john.lane@ksc.nasa.gov 

Subject: HAS Data Request: HAS001184684 Completed 

Your request (HAS0011 84684 ) has completed. The procedure to ftp the data is as 
follows : 

http: //wwwl . ncdc . noaa . qov/pub/has/HASOOl 18468 

NOTE: You must pick up your data within 7 days of this notice. 

Thank you, 

National Climatic Data Center 
E-mail: NCDC.Orders@noaa.gov 
Internet: www . ncdc . noaa . gov 

RADFILES.TXT 

The *.txt files downloaded from the NCDC ftp site are listed by the user, using a text editor, in 

RadFiles.txt. 
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1 . 

Pad-A 8x8k S.bmp 




2. 

28.5723 -80. 

6452 




3. 

7 000KMLB 

SDUS52 

N0RMLB 200702262214.txt 

10000.0 



4. 

7000KMLB 

SDUS22 

N1RMLB 200702262214. txt 

10000.0 



5. 

7000KMLB 

SDUS22 

N2RMLB 200702262214.txt 

10000.0 



6. 

7000KMLB 

S DUS32 

N3RMLB 200702262214 . txt 

10000.0 



7. 

# 






8. 

Zxyl.dat 


Zxy2 . dat 

40 

200 

10 

9. 

Zxzl .dat 


Zxz2 . dat 

40 

100 

4000 

10. 

Zyzl.dat 


Zyz2 . dat 

40 

100 

4000 


Using the above file as an example, the following notes can be generalized to other cases: 

■ Line 1 is a bitmap image supplied by the user of the x-y plane background, sized to a 200 x 

200 pixel image. The image represents a square plot region whose dimensions correspond 
to the product of the third and fourth entries in line 8 in meters. 

■ The lat/lon of the lower left comer of the plot region (southwest comer) is given in line 2. 

■ Lines 3 - 6 are the files downloaded via the NCDC ftp server described in the previous 

section. Note that not all four elevations are required. Any combination can be listed, but 
must be in vertical ascending order. The numerical parameter which follows the filename is 
the distance in meters from the region origin (given by line 2) that will define the region to 
be processed for plotting. Note that this value may be equal to or greater than the plot 
region spatial extent. 

■ Line 7 marks the end of the radar data file name entries. One to four entries are valid. 

■ Line 8 is the x-y plane parameters. The first string defines an output file that contains the x 

(north-south), y (east-west), and dB Z data, where dBZ is the radar reflectivity data. A spatial 
interpolation has been performed on the data contained in this file. The second file name 
corresponds to similar plot data, but without spatial interpolation. The third entry is x-y plot 
cell size in meters. The fourth entry is the number of cells, as discussed previously, the 
product of these two values defines the extent of the plot region in the x-y plane. The last 
parameter is the z position (elevation) in meters of the x-y plot plane. 

■ Line 9 is the x-z plane parameters. The first string defines an output file that contains the x 

(north-south), z (elevation), and dBZ data. A spatial interpolation has been performed on the 
data contained in this file. The second file name corresponds to similar plot data, but 
without spatial interpolation. The third entry is x-z plot cell size in meters. The fourth entry 
is the number of cells, where the product of these two values defines the extent of the plot 
region in the x-z plane. The last parameter is the y position (north-south position) in meters 
of the x-z plot plane. 

■ Line 10 is the y-z plane parameters. The first string defines an output file that contains they 

(east-west), z (elevation), and dBZ data. A spatial interpolation has been performed on the 
data contained in this file. The second file name corresponds to similar plot data, but 
without spatial interpolation. The third entry is y-z plot cell size in meters. The fourth entry 
is the number of cells, where the product of these two values defines the extent of the plot 
region in the y-z plane. The last parameter is the x position (east-west position) in meters of 
the y-z plot plane. 
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SITE COORD.DAT 


This input file contains the lat/lon of three disdrometer sites. There is no processing done in 
3DRadPlot requiring this data, other than small marks inserted in the 2D data plot files, 
corresponding to these site locations. 

Note: This file must exist with three entries, but the data contained within is not used by 
3DRadPlot. This file is not at all necessary or relevant unless disdrometer data is to be plotted 
by another application (not included here). 


OUTPUT FILES 

■ Arga.dat - this file corresponds exactly to the coordinates in Site Coord.dat but is represented 

by x-y values, in meters, from the origin of the x-y plot region, defined by line 1 of 
RadFiles.txt. 

■ Log.dat - this file contains miscellaneous information extracted from the Level III files, as 

well as processing comments. This file can be used for debugging. 

■ Six 2D plot output files, corresponding to the six names given in lines 8- 10 of RadFiles.txt. 

■ Pnnnn.dat - contains the x ,y, z, and t volume scan cell coordinates defined by the bounding 

box formed from the plot origin, line 2, and the region extent entry in lines 3-6 of 
RadFiles.txt. This file contains the 4D coordinates which would be used in the disdrometer 
interpolation/extrapolation for comparison to the radar data. 

■ Z3D.bmp - is a bitmap 3D plot of the 2D data. The colors represent dBZ increments of 5, 

starting with 5, ending with 65. 


EXECUTABLE 
The executable is: 

3DRadPlot.exe 

If all input files are in place, execution of the above file will generate the output files described 
above. 
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APPENDIX C 


4 th ORDER RUNGE-KUTTA FOR 2 nd ORDER DIFFERENTIAL EQUATIONS 


a, - v„_,A / 

Pi = X n-I + - 2 ®1 

y, =a(«A/,x I ,_ l ,v ( ,_ 1 )A/ 

8 i = v „-.+tY. 

a 2 =8, A / 

P 2 = X n-1 + 2 a 2 

1 2 =a((« + i)A/,p„8 l )A/ 

&2 =V „-I + TY2 

«3 = $2 A ' 

P 3 = X »-l + ®3 
Y3=a((n+i)A/,P 2 ,8 2 )A/ 

83 = V b-1+Y3 

«4=M' 

Y 4 =a((« + 1)A/,P 3 ,8 3 )A/ 

X » = X w— 1 + i(«| + 2 «2 + 2 «3 + « 4 ) 

v » = V B-i+i(Yi+2Y 2 + 2 Y3+Y 4 ) 

http://www.phvsics.orst.edu/~rubin/nacphv/ComPhvs/DIFFEO/mvdif2/ 
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